#include <TROOT.h>

void Z4F_p2FIT( int iOa2) 
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

	  //////////////////
	  // Input   ///////
	  //////////////////
	
	char Oa2print[8]; whichOa2(iOa2, Oa2print);
 
        char dir[1024];
	char rgiprint[1024];
	float theta[4];
	float cut;

	float **p2_interval=(float**)malloc(sizeof(float*)*(2)); // M1 & M2
	for (int im=0; im<2; im++) p2_interval[im]=(float*)malloc(sizeof(float)*(2));  // begin & end
	
	 // basic input file
	char input_dat[1024];
	sprintf(input_dat,"../INPUT/input_fit.dat");
	FILE *inputdat=fopen(input_dat,"r");
	fscanf(inputdat,"%i",&L);
        fscanf(inputdat,"%i",&T);
	fscanf(inputdat,"%f %f",&p2_interval[0][0], &p2_interval[0][1]);
        fscanf(inputdat,"%f %f",&p2_interval[1][0], &p2_interval[1][1]);
        fscanf(inputdat,"%f",&cut);
        fscanf(inputdat,"%i",&Nclust);
        fscanf(inputdat,"%s",dir);
        fscanf(inputdat,"%s",rgiprint);

	

   /////////////////////////
  //// Momenta set
  //////////////////////////////
  char momentafile[1024]; sprintf(momentafile,"../INPUT/mom.dat");
  read_momenta_data(momentafile);
  Equivalent_Momenta(SinP2);

 //////////////////////////////////
 // IO files
  ///////////////////////////////

  char file[1024];

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

  //////////////////////////////////
  // Read Zs
   /////////////////////////////

     double ****ZpPC=malloc(sizeof(double**)*(5));
     for (int iop1=0; iop1<5; iop1++) ZpPC[iop1]=malloc(sizeof(double**)*(5));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC[iop1][iop2]=malloc(sizeof(double*)*(Nmom_eq));
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) for (int imom=0; imom<Nmom_eq; imom++) ZpPC[iop1][iop2][imom]=malloc(sizeof(double)*(Nclust+1));

     double **Zq=malloc(sizeof(double*)*(Nmom_eq));
     for (int imom=0; imom<Nmom_eq; imom++) Zq[imom]=malloc(sizeof(double)*(Nclust+1));

     double *aux_clust=malloc(sizeof(double)*(Nclust+1));

     sprintf(file,"Zq%s_%s",rgiprint,Oa2print);

     for (int imom=0; imom<Nmom_eq; imom++)
     {
     for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
     {
                sprintf(file,"ZpPC%d%d%s_%s",iop1+1,iop2+1,rgiprint,Oa2print);
                read_input_onemomenta(dir,file,ZpPC[iop1][iop2][imom],imom);
     }

     read_input_onemomenta(dir,file,Zq[imom],imom);
     }
   //////////////////////////////
   // p2 fit
   ////////////////////////////////
	
  // auxiliar for minuit routines
  // in practice I'll need Nmom_for_fit < Nmom 
 double *ZAUX=malloc(sizeof(double)*(Nmom));
  double **ZAUX_CLUST=malloc(sizeof(double*)*(Nmom));
  for (int imom=0; imom<Nmom; imom++) ZAUX_CLUST[imom]=malloc(sizeof(double)*(Nclust+1));

  double *P2AUX=malloc(sizeof(double)*(Nmom));

  double ****ZpPC_FIT=malloc(sizeof(double***)*(2));
  for (int im=0; im<2; im++) ZpPC_FIT[im]=malloc(sizeof(double**)*(5));
  for (int im=0; im<2; im++) for (int iop1=0; iop1<5; iop1++) ZpPC_FIT[im][iop1]=malloc(sizeof(double*)*(5));
  for (int im=0; im<2; im++) for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC_FIT[im][iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

  double **Zq_FIT=malloc(sizeof(double*)*(2));
  for (int im=0; im<2; im++) Zq_FIT[im]=malloc(sizeof(double)*(Nclust+1));


  double ****ZpPC_SLOPE=malloc(sizeof(double***)*(2));
  for (int im=0; im<2; im++) ZpPC_SLOPE[im]=malloc(sizeof(double**)*(5));
  for (int im=0; im<2; im++) for (int iop1=0; iop1<5; iop1++) ZpPC_SLOPE[im][iop1]=malloc(sizeof(double*)*(5));
  for (int im=0; im<2; im++) for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC_SLOPE[im][iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

  double **Zq_SLOPE=malloc(sizeof(double*)*(2));
  for (int im=0; im<2; im++) Zq_SLOPE[im]=malloc(sizeof(double)*(Nclust+1));

  double ****ZpPC_CHI2=malloc(sizeof(double***)*(2));
  for (int im=0; im<2; im++) ZpPC_CHI2[im]=malloc(sizeof(double**)*(5));
  for (int im=0; im<2; im++) for (int iop1=0; iop1<5; iop1++) ZpPC_CHI2[im][iop1]=malloc(sizeof(double*)*(5));
  for (int im=0; im<2; im++) for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++) ZpPC_CHI2[im][iop1][iop2]=malloc(sizeof(double)*(Nclust+1));

  double **Zq_CHI2=malloc(sizeof(double*)*(2));
  for (int im=0; im<2; im++) Zq_CHI2[im]=malloc(sizeof(double)*(Nclust+1));



  int Nmom_for_fit;	

	// zq
   for (int im=0; im<2; im++)
      for (int iclust=0; iclust<Nclust+1; iclust++)
      {
          int imom_for_fit=0;
          for (int imom=0; imom<Nmom_eq; imom++)
          {
                if (SinP2_eq[imom] > p2_interval[im][0]  && SinP2_eq[imom]< p2_interval[im][1] && SinP4_eq[imom]/pow(SinP2_eq[imom],2)<cut )
                {
                         ZAUX[imom_for_fit]=Zq[imom][iclust];
                         for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imom_for_fit][iclust_loop]=Zq[imom][iclust_loop];
                         P2AUX[imom_for_fit]=SinP2_eq[imom];
                         imom_for_fit++;
                   }
          }
          Nmom_for_fit=imom_for_fit;
          // M1
          if (im==0) lin_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Zq_FIT[im][iclust],Zq_SLOPE[im][iclust],Zq_CHI2[im][iclust]); 
           // M2
           if (im==1) {wprom_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,Zq_FIT[im][iclust],Zq_CHI2[im][iclust]); Zq_SLOPE[im][iclust]=0.;}

      }

	// Z4f
   for (int im=0; im<2; im++)
      for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
      for (int iclust=0; iclust<Nclust+1; iclust++)
      {
          int imom_for_fit=0;
          for (int imom=0; imom<Nmom_eq; imom++)
          {
                if (SinP2_eq[imom] > p2_interval[im][0]  && SinP2_eq[imom]< p2_interval[im][1] && SinP4_eq[imom]/pow(SinP2_eq[imom],2)<cut )
                {
                         ZAUX[imom_for_fit]=ZpPC[iop1][iop2][imom][iclust];
                         for (int iclust_loop=0; iclust_loop<Nclust+1; iclust_loop++) ZAUX_CLUST[imom_for_fit][iclust_loop]=ZpPC[iop1][iop2][imom][iclust_loop];
                         P2AUX[imom_for_fit]=SinP2_eq[imom];
                         imom_for_fit++;
                   }
          }
          Nmom_for_fit=imom_for_fit;

          // M1
          if (im==0) lin_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,ZpPC_FIT[im][iop1][iop2][iclust],ZpPC_SLOPE[im][iop1][iop2][iclust],ZpPC_CHI2[im][iop1][iop2][iclust]);
           // M2
          if (im==1) {wprom_fit(ZAUX,ZAUX_CLUST,P2AUX,Nmom_for_fit,Nclust,ZpPC_FIT[im][iop1][iop2][iclust],ZpPC_CHI2[im][iop1][iop2][iclust]); ZpPC_SLOPE[im][iop1][iop2][iclust]=0.;}
      }


	// write output

	double **aux_allmomenta=malloc(sizeof(double*)*(Nmom));
	for (int imom=0; imom<Nmom; imom++) aux_allmomenta[imom]=malloc(sizeof(double)*(Nclust+1));

	for (int im=0; im<2; im++)
	{
		sprintf(file,"%s/p2FIT/M%d/Zq_%s%s",dir,im+1,Oa2print,rgiprint);
		write_output_fitmomenta_real(file,Zq_FIT[im]);
                sprintf(file,"%s/p2FIT/M%d/Zq_%s%s_slope",dir,im+1,Oa2print,rgiprint);
                write_output_fitmomenta_real(file,Zq_SLOPE[im]);
		
		for (int iop1=0; iop1<5; iop1++) for (int iop2=0; iop2<5; iop2++)
		{
			sprintf(file,"%s/p2FIT/M%d/ZpPC%d%d_%s%s",dir,im+1,iop1+1,iop2+1,Oa2print,rgiprint);
			write_output_fitmomenta_real(file,ZpPC_FIT[im][iop1][iop2]);
                        sprintf(file,"%s/p2FIT/M%d/ZpPC%d%d_%s%s_slope",dir,im+1,iop1+1,iop2+1,Oa2print,rgiprint);
                        write_output_fitmomenta_real(file,ZpPC_SLOPE[im][iop1][iop2]);
		}
			
	}

}
